﻿#define _CRT_SECURE_NO_WARNINGS

#include <iostream>
#include <string>
#include "gdal_priv.h"
#include <cassert>
#include "boost/filesystem.hpp"
#include <sstream>

using namespace std;

template<typename T>
void tiles(GDALDataset* ds, GDALDataType dataType, double filldata, string outdir);
const int TILE_SIZE = 256;
const int TILE_SIZE_2 = TILE_SIZE * TILE_SIZE;
const int MAX_LEVEL = 18;
const double RESO_ARR[MAX_LEVEL + 1] = { 1.4062500,0.7031250,0.3515625,0.1757813,0.0878906,
                            0.0439453,0.0219727,0.0109863,0.0054932,0.0027466,
                            0.0013733,0.0006866,0.0003433,0.0001717,0.0000858,
                            0.0000429,0.0000215,0.0000107,0.0000054 };


//level from 0 to MAX_LEVEL.
int calculateLevel(double inputResolution) {
    
    if (inputResolution <= 0) return 0;
    if (inputResolution >= RESO_ARR[0]) return 0;
    if (inputResolution <= RESO_ARR[MAX_LEVEL]) return MAX_LEVEL;
    for (int i = 0; i < MAX_LEVEL; ++i)
    {
        if (inputResolution < RESO_ARR[i] && inputResolution >= RESO_ARR[i + 1]) {
            double mid = (RESO_ARR[i + 1]+ RESO_ARR[i]) * 0.5;
            if (inputResolution >= mid) return i;
            else return i + 1;
        }
    }
    return -1; // this is unreachable!
}

int main(int argc, char* argv[])
{
    cout << "This is a very simple program to cut tiles for WGS84 geotiff datas. Each output tile is 256x256 BSQ binary data. In level 0, just one tile and the upper half is valid." << endl;
    cout << "This program is a very simple utility of OEE for cutting demo data into tiles. 2024/4/4/" << endl;
    cout << "Input tif file should be WGS84 spatial reference. Support max level is 18 ( about 0.59545m/px 0.000005364418 in degree )" << endl;
    cout << "outdir and outtiles are /{outdir}/{z}/{x}/{y}.raw" << endl;
    cout << "version: v1.0.0.1" << endl;
    cout << "usage: oee_tiles /some/dir/file.tif filldata /out/put/dir" << endl;
    GDALAllRegister();

    if (argc != 4) {
        cout << "argc not 4, argc=" << argc << endl;
        return 11;
    }
    string infile = argv[1];
    double filldata = atof(argv[2]);
    string outdir = argv[3];

    cout << "infile " << infile << endl;
    cout << "filldata " << filldata << endl;
    cout << "outdir " << outdir << endl;

    bool out1ok = boost::filesystem::exists(outdir.c_str());
    if (out1ok == false) {
        cout << outdir << " not exists, try to make it..." << endl;
        bool mkdir1ok = boost::filesystem::create_directories(outdir.c_str());
        if (mkdir1ok) {
            cout << "mkdir ok" << endl; 
        }
        else {
            cout << "Failed to mkdir of " << outdir << endl;
            return 12;
        }
    }

    double trans[6];
    GDALDataset* ds = (GDALDataset*)GDALOpen(infile.c_str(), GA_ReadOnly);
    if (ds == 0) {
        cout << "Failed to open input file " << infile << endl; 
        return 13;
    }
    const int xsize = ds->GetRasterXSize();
    const int ysize = ds->GetRasterYSize();
    const int nbands = ds->GetRasterCount();
    cout << "raster xsize ysize nbands " << xsize << " " << ysize << " " << nbands << endl;
    ds->GetGeoTransform(trans);
    GDALDataType dataType = ds->GetRasterBand(1)->GetRasterDataType();
    if (dataType == GDT_Byte) {      
        tiles<unsigned char>(ds, GDT_Byte,filldata,outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_UInt16) {
        tiles<unsigned short>(ds, GDT_UInt16, filldata, outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_Int16) {
        tiles<short>(ds, GDT_Int16, filldata, outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_UInt32) {
        tiles<unsigned int>(ds, GDT_UInt32, filldata, outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_Int32) {
        tiles<int>(ds, GDT_Int32, filldata, outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_Float32) {
        tiles<float>(ds, GDT_Float32, filldata, outdir);
        GDALClose(ds);
    }
    else if (dataType == GDT_Float64) {
        tiles<double>(ds, GDT_Float64, filldata, outdir);
        GDALClose(ds);
    }
    else {
        cout << "unsupported data type " << dataType << endl;
        GDALClose(ds);
        return 20;
    }
    cout << "done" << endl;
    return 0;
}

string int2str(int val) {
    stringstream ss;
    ss << val;
    return ss.str();
}

template<typename T>
void tiles(GDALDataset* ds,GDALDataType dataType, double filldata, string outdir)
{
    const int xsize = ds->GetRasterXSize();
    const int ysize = ds->GetRasterYSize();
    const int nbands = ds->GetRasterCount();
    double trans[6];
    ds->GetGeoTransform(trans);
    int theLevel = calculateLevel(trans[1]);
    assert(theLevel >= 0);
    cout << "theLevel " << theLevel << endl;
    double left = trans[0];
    double right = trans[0] + xsize * trans[1];
    double top = trans[3];
    double bottom = trans[3] + ysize * trans[5];

    vector< vector<T> > bandDataArray(nbands);
    vector<T*> bandDataPtrArray(nbands);
    for (int ib = 0; ib < nbands; ++ib) {
        bandDataArray[ib].resize(xsize * ysize);
        ds->GetRasterBand(ib + 1)->RasterIO(GF_Read,
            0, 0, xsize, ysize, bandDataArray[ib].data(),
            xsize, ysize, dataType, 0, 0, 0);
        bandDataPtrArray[ib] = bandDataArray[ib].data();
    }
    boost::filesystem::path rootdirpath(outdir.c_str());
    vector<T> tt(TILE_SIZE * TILE_SIZE*nbands);
    for (int ilevel = 0; ilevel <= theLevel; ++ilevel)
    {
        boost::filesystem::path zdir = rootdirpath / int2str(ilevel);
        bool mkdirz = boost::filesystem::create_directories(zdir);
        if (mkdirz == false) {
            cout << "failed to mkdir for z " << zdir << endl;
            exit(21);
        }
        int maxxtilecount = pow(2, ilevel);
        int maxytilecount = maxxtilecount / 2;
        if (maxytilecount == 0) maxytilecount = 1;
        maxxtilecount -= 1;
        maxytilecount -= 1; 

        int tilex0 = ((int)( (left+180.0) / RESO_ARR[ilevel])) / TILE_SIZE;
        int tilex1 = MIN(maxxtilecount, ceil( (int)( (right+180.0 - RESO_ARR[ilevel]/2) / RESO_ARR[ilevel]) / TILE_SIZE)) ;

        int tiley0 = ((int)( (90.0-top) / RESO_ARR[ilevel])) / TILE_SIZE;
        int tiley1 = MIN(maxytilecount,ceil((int)( (90.0-bottom) / RESO_ARR[ilevel]) / TILE_SIZE));

        double lreso = RESO_ARR[ilevel];
        double tile_dReso = TILE_SIZE * lreso;
        double hreso = lreso * 0.5; //half it

        for (int tilex = tilex0; tilex <=  tilex1; ++tilex) {
            boost::filesystem::path xdir = zdir / int2str(tilex);
            double tilelonUL = -180 + tile_dReso * tilex + hreso;//center leftup

            bool mkdirx = boost::filesystem::create_directories(xdir);
            if (mkdirx == false) {
                cout << "failed to mkdir for x " << xdir << endl;
                exit(22);
            }

            for (int tiley = tiley0; tiley <= tiley1; ++tiley) {
                double tilelatUL = 90 - tile_dReso * tiley - hreso;//center leftup 
                boost::filesystem::path tilefilename = xdir / (int2str(tiley)+".raw");
                std::fill(tt.begin(), tt.end(), (T)filldata);
                double pixel_y = tilelatUL;
                T* ttptr = tt.data();
                for (int iy = 0; iy < TILE_SIZE; ++iy) {
                    double pixel_x = tilelonUL;
                    int input_y = (pixel_y - trans[3]) / trans[5];
                    if (input_y >= 0 && input_y < ysize) {
                        for (int ix = 0; ix < TILE_SIZE; ++ix) {
                            int input_x = (pixel_x - trans[0]) / trans[1];
                            if (input_x >= 0 && input_x < xsize) {
                                for (int ib = 0; ib < nbands; ++ib)
                                    *(ttptr + TILE_SIZE_2 * ib + iy * TILE_SIZE + ix) = *(bandDataPtrArray[ib] + input_y * xsize + input_x);
                            }
                            pixel_x += lreso;
                        }
                    }
                    pixel_y -= lreso;
                }
                //write down files.
                {
                    FILE* pftile = fopen(tilefilename.string().c_str(), "wb");
                    if (pftile == 0) {
                        cout << "failed to create tile file " << tilefilename << endl;
                        exit(23);
                    }
                    fwrite(tt.data(), 1, TILE_SIZE_2 * nbands * sizeof(T), pftile);
                    fclose(pftile);
                }
            }
        }
    }
}